/**
 * @file Multigrid_Solver.h
 * @brief 
 * @author Linhuo (1066117119@qq.com)
 * @version 1.0
 * @date 2021-06-08
 * 
 * @copyright Copyright (c) 2021  Linhuo
 */
#ifndef _PRO2_MULTIGRID_SOLVER_H_
#define _PRO2_MULTIGRID_SOLVER_H_
#include <iostream>
#include <eigen3/Eigen/Core>
#include <functional>
#include <cmath>
#include "Intp_Solver.h"
#include "Intp_Finder.h"
#include "Restric_Solver.h"
#include "Restrict_Finder.h"
typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> MatrixXd;
typedef Eigen::Matrix<double, Eigen::Dynamic, 1> VectorXd;
typedef std::function<double(double *)> Function;

/**
 * @brief An enumeration of Smoothing Function；
 * 		  WeightedJacobi ,omega = 2/3　= 　1;
 * 		　GaussSeidel 　= 　2；
*/
enum Smoother
{
    WeightedJacobi = 1,
    GaussSeidel = 2
};

/**
 * @brief A pure virtual base class with primary connector and data structure;
 * A base for all concrete implementation for Multigrid Solver.
 */
class Multigrid_Solver
{
protected:
    /**
     * @brief The dimension of PDE.
     */
    int _dimension = 1;
    /**
     * @brief The finest grid number (for one side). (Preference: 128)
     */
    int _num_of_grid = 128;
    /**
     * @brief  Store original iterative initial value. ( with ghost cell)
     */
    VectorXd _u_origin;
    /**
     * @brief Store original right-hand-sentence. ( with ghost cell)
     */
    VectorXd _f_origin;
    /**
     * @brief Store boundary condition.
     */
    MatrixXd _BC;
    /**
     * @brief Store smoothing type. (Preference: GS method)
     */
    Smoother SmoothType = GaussSeidel;
    /**
     * @brief Store restriction type. (Preference: Full-Weighting)
     */
    Logo_Restrict RestricType = FullWeighting;
    /**
     * @brief Store interpolation type. (Preference: Linear)
     */
    Logo_Intper IntpType = Linear;

public:
    /**
     * @brief Get the f_origin object.
     * @return VectorXd 
     */
    const VectorXd get_f_origin() { return _f_origin; };
    /**
     * @brief Get the u_origin object.
     * @return VectorXd 
     */
    const VectorXd get_u_origin() { return _u_origin; };
    /**
     * @brief Get the Boundrary Condition.
     * @return MatrixXd 
     */
    const MatrixXd get_BC() { return _BC; };
    /**
     * @brief Get the h object.
     * @return double 
     */
    const double get_h() { return 1.0 / _num_of_grid; };
    /**
     * @brief Get the dimension object.
     * @return const int 
     */
    const int get_dimension(){return _dimension;};
    /**
     * @brief Set the Smooth Type object.
     * @param  SmoothTypeCode 
     */
    void setSmoothType(int SmoothTypeCode) { SmoothType = (enum Smoother)SmoothTypeCode; };
    /**
     * @brief Set the Restirc Type object.
     *  FullWeighting 　= 　1； Injection 　= 　2；
     * @param  RestricTypeCode 
     */
    void setRestircType(int RestricTypeCode) { RestricType = (enum Logo_Restrict)RestricTypeCode; };
    /**
     * @brief Set the Intp Type object.
     *  Linear 　= 　1；Quadratic 　= 　2；
     * @param  IntpTypeCode  
     */
    void setIntpType(int IntpTypeCode) { IntpType = (enum Logo_Intper)IntpTypeCode; };
    /**
     * @brief Set the Num Of Grid object.
     * @param  num          
     */
    void setNumOfGrid(int num) { _num_of_grid = num; };
    /**
     * @brief Set the Dimension object
     * @param  dim              
     */
    void setDimension(int dim) {_dimension = dim;};
    /**
     * @brief Initialize numofgrid, initial value, and boundary condition.(with ghost cell)
     * @param  num             numofgrid
     * @param  _F              Rhs Function.
     * @param  _BCF             Boundary Condition Function.
     */
    void Initialize(int num, Function _F, Function _BCF);
    /**
     * @brief Concretely initialize initial value.(with ghost cell)
     */
    virtual void Initialize_u() = 0;
    /**
     * @brief  Concretely initialize boundary condition.(with ghost cell)
     * @param  _BCF             Boundary Condition Function
     */
    virtual void Initialize_BC(Function _BCF) = 0;
    /**
     * @brief Concretely initialize RHS.(with ghost cell)
     * @param  _F               RHS Function
     */
    virtual void InitializeRhs(Function _F) = 0;
    /**
     * @brief Specific Smoothing Solver for k-Dimension Possion Equation with Dirichlet Boundary.
     * @param  _u               Initial Value (with ghost cell)
     * @param  _f               Boundary Condition (w.g.c.)
     * @param  _h               Width of Finest Grid (w.g.c.)
     * @return VectorXd Smoothed Result
     */
    virtual VectorXd SmoothingforVC(VectorXd _u, VectorXd _f, double _h) = 0;
    /**
     * @brief Specific Residual Solver for k-Dimension Possion Equation with Dirichlet Boundary.
     * @param  _u               Initial Value (with ghost cell)
     * @param  _f               Boundary Condition (w.g.c.)
     * @param  _h               Width of Up-dated Grid (w.g.c.)
     * @return VectorXd Residual
     */
    virtual VectorXd ResidualSolver(const VectorXd _u, const VectorXd _f, const double _h) = 0;
    /**
     * @brief Specific Restrictor for k-Dimension vector.
     * @param  _r               
     * @return VectorXd Restricted Vector
     */
    // virtual VectorXd RestrictSolver(VectorXd _r) = 0;
    VectorXd RestrictSolver(VectorXd _r);
    /**
     * @brief Specific Interpolater for k-Dimension vector.
     * @param  _r               
     * @return VectorXd Interpolated Vector
     */
    VectorXd IntpSolver(VectorXd _r);
    /**
     * @brief A Single Vcycle-Solver, i.e. with (\mu_1 = 1, \mu_2 = 1).
     * @param  _u               Initial Value (with ghost cell)
     * @param  _f               Boundary Condition (w.g.c.)
     * @param  _h               Width of Up-dated Grid (w.g.c.)
     * @return VectorXd Solution
     */
    VectorXd Vcycle(VectorXd _u, VectorXd _f, double _h);
    /**
     * @brief Specific Preparation for FMG.
     * @param  _f               Boundary Condition.
     * @return VectorXd Initial Value.(with ghost cell)
     */
    virtual VectorXd PrepareFMG(VectorXd _f) = 0;
    /**
     * @brief Get the now num grid object
     * @param  _f               Vector
     * @return int now num gird
     */
    virtual int get_now_num_grid(VectorXd _f) = 0;
    /**
     * @brief FMG cycle with single V-cycles.
     * @param  _f               Boundary Condition
     * @param  _h               Width of Up-dated Grid
     * @return VectorXd 
     */
    VectorXd FMGcycle(VectorXd _f, double _h);
};

void Multigrid_Solver::Initialize(int num, Function _F, Function _BCF)
{
    setNumOfGrid(num);
    Initialize_BC(_BCF);
    InitializeRhs(_F);
    Initialize_u();
};

VectorXd Multigrid_Solver::Vcycle(VectorXd _u, VectorXd _f, double _h)
{
    VectorXd __u, __r, __eps;
    __u = SmoothingforVC(_u, _f, _h);
    __r = ResidualSolver(__u, _f, _h);
    __r = RestrictSolver(__r);
    __eps.setZero(__r.rows());
    if (get_now_num_grid(__eps) == 2)
    {
        __eps = SmoothingforVC(__eps, __r, 2 * _h);
    }
    else
    {
        __eps = Vcycle(__eps, __r, 2 * _h);
    };
    __u = __u + IntpSolver(__eps);
    __u = SmoothingforVC(__u, _f, _h);
    return __u;
};

VectorXd Multigrid_Solver::FMGcycle(VectorXd _f, double _h)
{
    VectorXd __u, __r;
    int N;
    N = get_now_num_grid(_f);
    __u = PrepareFMG(_f);
    if (N == 2)
    {
        __u = SmoothingforVC(__u, _f, _h);
    }
    else
    {
        __r = RestrictSolver(_f);
        __u = FMGcycle(__r, 2 * _h);
        __u = IntpSolver(__u);
        __u = Vcycle(__u, _f, _h);
    };
    return __u;
};

VectorXd Multigrid_Solver::IntpSolver(VectorXd _r)
{
    Intp_Finder helper(IntpType,_dimension);
    return helper.Interpolate(_r);
};

VectorXd Multigrid_Solver::RestrictSolver(VectorXd _r)
{  
    Restrict_Finder helper(RestricType,_dimension);
    return helper.Restrict(_r);
};

#else
//DO NOTHING
#endif